from astropy.io import fits
from astropy.table import Table
import numpy as np
import matplotlib.pyplot as plt
from astropy.stats import sigma_clipped_stats
from photutils import datasets
from photutils import CircularAperture
from photutils import DAOStarFinder
hdu = datasets.load_star_image('quadRU.fits')
data = hdu.data#[0:401, 0:401]
mean, median, std = sigma_clipped_stats(data, sigma=3.0)
print((mean, median, std))
(3770.4896097401975, 3744.0, 237.21215226532127)
#Plotting the image
fig, ax = plt.subplots(1,2)
fig.set_size_inches(16,12)
ax[0].imshow(data,cmap=plt.cm.gray) #change colour to gray => cmap=plt.cm.gray
ax[1].imshow(data[401:801, 401:801] ,cmap=plt.cm.gray)
print("The number of pixels is {0}".format(data.size))
print("The mean is {0} and the STD ".format(np.mean(data)),format(np.std(data)))
The number of pixels is 1121481 The mean is 3988.5913154123878 and the STD 1069.1249436838032
#STEP 6
daofind = DAOStarFinder(fwhm=3.0, threshold=27.7*std) #Full width at half maximum(fwhm) // good value = 27.7
sources = daofind(data)# - median)
for col in sources.colnames:
sources[col].info.format = '%.8g' # for consistent table output
sources.sort('xcentroid')
sources
| id | xcentroid | ycentroid | sharpness | roundness1 | roundness2 | npix | sky | peak | flux | mag |
|---|---|---|---|---|---|---|---|---|---|---|
| int64 | float64 | float64 | float64 | float64 | float64 | int64 | float64 | float64 | float64 | float64 |
| 170 | 0.91894685 | 840.89515 | 0.30017139 | 0.42374031 | 0.73439404 | 25 | 0 | 12645 | 1.4511468 | -0.40427836 |
| 132 | 1.3110953 | 603.52503 | 0.3070012 | -0.17534705 | 0.8687376 | 25 | 0 | 12586 | 1.0172121 | -0.018528787 |
| 48 | 1.3819622 | 194.14896 | 0.34887691 | -0.040359472 | 0.91937352 | 25 | 0 | 12439 | 1.0000395 | -4.2841351e-05 |
| 194 | 13.709729 | 966.33218 | 0.37653395 | -0.061260617 | -0.0047948456 | 25 | 0 | 12551 | 1.0262129 | -0.028093697 |
| 28 | 14.85393 | 104.58858 | 0.41325279 | -0.15754959 | 0.0014790695 | 25 | 0 | 12176 | 1.078323 | -0.081872155 |
| 143 | 15.308049 | 662.25654 | 0.40410494 | 0.092247523 | -0.21561658 | 25 | 0 | 12734 | 1.0472938 | -0.050171342 |
| 201 | 26.026207 | 990.82951 | 0.40973918 | -0.12748799 | -0.1263963 | 25 | 0 | 12289 | 1.1639493 | -0.16483515 |
| 41 | 31.916545 | 171.79805 | 0.40947179 | 0.073722424 | -0.097077158 | 25 | 0 | 12696 | 1.1165765 | -0.11972121 |
| 63 | 35.353299 | 267.10879 | 0.40362557 | -0.1011558 | -0.10080115 | 25 | 0 | 12498 | 1.0806916 | -0.084254463 |
| 21 | 37.115099 | 83.290352 | 0.55358402 | 0.018959982 | 0.032717845 | 25 | 0 | 11530 | 1.007856 | -0.00849618 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 151 | 974.21449 | 695.40998 | 0.48341372 | -0.31325491 | 0.011266162 | 25 | 0 | 11935 | 1.0067901 | -0.0073473202 |
| 5 | 981.81682 | 32.809747 | 0.49815515 | -0.038007005 | -0.15331728 | 25 | 0 | 11682 | 1.0966799 | -0.10019975 |
| 111 | 1014.2458 | 456.11819 | 0.38251006 | 0.075475958 | 0.027730015 | 25 | 0 | 12529 | 1.0789151 | -0.082468222 |
| 125 | 1023.7488 | 564.6262 | 0.44042157 | 0.30361012 | -0.042326908 | 25 | 0 | 11600 | 1.0066716 | -0.007219594 |
| 114 | 1039.4274 | 492.84468 | 0.43837681 | 0.031428091 | -0.24080226 | 25 | 0 | 12320 | 1.0214756 | -0.023070039 |
| 27 | 1041.6488 | 101.70837 | 0.45402399 | 0.19568416 | -0.058768625 | 25 | 0 | 11825 | 1.0237259 | -0.025459175 |
| 35 | 1043.7634 | 119.08471 | 0.49411584 | 0.031467285 | -0.10235053 | 25 | 0 | 12143 | 1.1573171 | -0.15863093 |
| 138 | 1044.666 | 632.90166 | 0.53714543 | -0.1248506 | -0.12651939 | 25 | 0 | 11895 | 1.0521495 | -0.055193654 |
| 167 | 1044.7256 | 804.82416 | 0.53046292 | -0.11265257 | -0.044356247 | 25 | 0 | 11594 | 1.0156725 | -0.016884211 |
| 96 | 1048.2053 | 380.56477 | 0.41206942 | -0.19084793 | -0.080891852 | 25 | 0 | 12445 | 1.0497809 | -0.052746648 |
import matplotlib.pyplot as plt
from astropy.visualization import SqrtStretch
from astropy.visualization.mpl_normalize import ImageNormalize
from photutils import CircularAperture
plt.figure(figsize=(20,20))
positions = (sources['xcentroid'], sources['ycentroid'])
apertures = CircularAperture(positions, r=4.)
norm = ImageNormalize(stretch=SqrtStretch())
plt.imshow(data, cmap='Greys', origin='lower', norm=norm)
apertures.plot(color='blue', lw=1.5, alpha=0.5)
#STEP 8
import matplotlib.pyplot as plt
from astropy.visualization import SqrtStretch
from astropy.visualization.mpl_normalize import ImageNormalize
from photutils import CircularAperture
from photutils import aperture_photometry
plt.figure(figsize=(20,20))
positions = (sources['xcentroid'], sources['ycentroid'])
apertures = CircularAperture(positions, r=4.)
phot_table = aperture_photometry(data, apertures, method='exact',subpixels=5)
norm = ImageNormalize(stretch=SqrtStretch())
plt.imshow(data, cmap='Greys', origin='lower', norm=norm)
apertures.plot(color='blue', lw=1.5, alpha=0.5)
phot_table = aperture_photometry(data, apertures, method='exact',subpixels=5)
phot_table['aperture_sum'].info.format = '%.8g' # for consistent table output
print(phot_table)
id xcenter ycenter aperture_sum
pix pix
--- ------------------ ------------------ ------------
1 0.9189468499910883 840.8951457248296 235087.24
2 1.311095259132758 603.5250306712844 387140.03
3 1.3819621776452984 194.1489634549567 372487.51
4 13.709729099697494 966.3321753368325 298805.15
5 14.853929564441659 104.5885806490168 275656.75
6 15.308049034659224 662.2565373858495 320886.39
7 26.026207222819796 990.8295092551749 271741.42
8 31.91654474923815 171.79804519605577 283941.14
9 35.353299330115014 267.1087864002644 285376.66
10 37.11509862938315 83.29035164102501 249927.49
... ... ... ...
199 974.2144947330132 695.4099754883719 265665.99
200 981.8168232536541 32.80974670705418 247157.03
201 1014.24581863924 456.1181928875807 293596.48
202 1023.7488485932108 564.6262045991909 266480.28
203 1039.4273600411575 492.84467803095777 287349.87
204 1041.648750654047 101.7083747902434 267495.91
205 1043.763408945062 119.08471405894595 254575.88
206 1044.666008741653 632.9016640068996 256337.64
207 1044.725609535972 804.8241634613814 253547.55
208 1048.2053425810136 380.56477025247375 292795.78
Length = 208 rows
phot_table = aperture_photometry(data, apertures, method='subpixel',subpixels=5)
phot_table['aperture_sum'].info.format = '%.8g' # for consistent table output
print(phot_table)
id xcenter ycenter aperture_sum
pix pix
--- ------------------ ------------------ ------------
1 0.9189468499910883 840.8951457248296 235553.8
2 1.311095259132758 603.5250306712844 388696.92
3 1.3819621776452984 194.1489634549567 372579.88
4 13.709729099697494 966.3321753368325 299362.88
5 14.853929564441659 104.5885806490168 275714.36
6 15.308049034659224 662.2565373858495 321257.12
7 26.026207222819796 990.8295092551749 271478.6
8 31.91654474923815 171.79804519605577 283271.96
9 35.353299330115014 267.1087864002644 285431.44
10 37.11509862938315 83.29035164102501 250929.24
... ... ... ...
199 974.2144947330132 695.4099754883719 265225.84
200 981.8168232536541 32.80974670705418 246734.24
201 1014.24581863924 456.1181928875807 293849.52
202 1023.7488485932108 564.6262045991909 266887.52
203 1039.4273600411575 492.84467803095777 287735.48
204 1041.648750654047 101.7083747902434 267436.12
205 1043.763408945062 119.08471405894595 254792.08
206 1044.666008741653 632.9016640068996 256692.12
207 1044.725609535972 804.8241634613814 253284.56
208 1048.2053425810136 380.56477025247375 292364.44
Length = 208 rows
phot_table = aperture_photometry(data, apertures, method='center',subpixels=5)
phot_table['aperture_sum'].info.format = '%.8g' # for consistent table output
print(phot_table)
id xcenter ycenter aperture_sum
pix pix
--- ------------------ ------------------ ------------
1 0.9189468499910883 840.8951457248296 226283
2 1.311095259132758 603.5250306712844 404626
3 1.3819621776452984 194.1489634549567 371062
4 13.709729099697494 966.3321753368325 298020
5 14.853929564441659 104.5885806490168 274731
6 15.308049034659224 662.2565373858495 326327
7 26.026207222819796 990.8295092551749 266936
8 31.91654474923815 171.79804519605577 283197
9 35.353299330115014 267.1087864002644 288381
10 37.11509862938315 83.29035164102501 252774
... ... ... ...
199 974.2144947330132 695.4099754883719 260311
200 981.8168232536541 32.80974670705418 254090
201 1014.24581863924 456.1181928875807 296865
202 1023.7488485932108 564.6262045991909 265634
203 1039.4273600411575 492.84467803095777 286708
204 1041.648750654047 101.7083747902434 266108
205 1043.763408945062 119.08471405894595 253822
206 1044.666008741653 632.9016640068996 259520
207 1044.725609535972 804.8241634613814 256679
208 1048.2053425810136 380.56477025247375 291884
Length = 208 rows
#STEP 10 - 11
from photutils import CircularAnnulus
apertures4 = CircularAperture(positions, r=3.)
annulus_apertures4 = CircularAnnulus(positions ,r_in=6 ,r_out=9)
apers4 = [apertures4, annulus_apertures4]
phot_table4 = aperture_photometry(data, apers4)
for col in phot_table4.colnames:
phot_table4[col].info.format = '%.8g' # for consistent table out
print(phot_table4)
#The aperture_sum_0 column refers to the first aperture (circular)
#the aperture_sum_1 column refers to the second aperture (annulus)
id xcenter ycenter aperture_sum_0 aperture_sum_1
pix pix
--- ---------- --------- -------------- --------------
1 0.91894685 840.89515 175857.43 299406.7
2 1.3110953 603.52503 283355.7 342294.06
3 1.3819622 194.14896 280215.01 317289.96
4 13.709729 966.33218 206998.37 524160.52
5 14.85393 104.58858 187510.09 533347.69
6 15.308049 662.25654 210624.24 594384.47
7 26.026207 990.82951 183965.79 521889.02
8 31.916545 171.79805 194901.82 518105.3
9 35.353299 267.10879 195810.43 514602.35
10 37.115099 83.290352 160962 535694.26
... ... ... ... ...
199 974.21449 695.40998 178183.94 537045.49
200 981.81682 32.809747 162077.25 522570.44
201 1014.2458 456.11819 200126.95 548544.89
202 1023.7488 564.6262 177732.56 544490.79
203 1039.4274 492.84468 194633.15 537354.49
204 1041.6488 101.70837 178287.04 524084.66
205 1043.7634 119.08471 168338.84 520448.17
206 1044.666 632.90166 166209.24 550289.95
207 1044.7256 804.82416 164841.21 536903.28
208 1048.2053 380.56477 198785.18 547217.71
Length = 208 rows
bkg_mean = phot_table4['aperture_sum_1'] / annulus_apertures.area()
bkg_sum = bkg_mean * apertures4.area()
final_sum = phot_table4['aperture_sum_0'] - bkg_sum
phot_table4['residual_aperture_sum'] = final_sum
phot_table4['residual_aperture_sum'].info.format = '%.8g' # for consistent table output
print(phot_table4['residual_aperture_sum'])
residual_aperture_sum
---------------------
154300.15
258710.53
257370.14
169258.81
149109.06
167828.56
146389.78
157598.24
158759.06
122392.01
...
120263.39
139516.67
124452.18
160631.72
138529.23
155943.63
140552.95
130866.57
126588.36
126184.18
159385.51
Length = 208 rows
from photutils import CircularAperture, CircularAnnulus
from astropy.visualization import simple_norm
norm = simple_norm(data, 'sqrt', percent=99)
plt.figure(figsize=(30,30))
plt.imshow(data, norm=norm)
apertures4.plot(color='white', lw=2)
annulus_apertures4.plot(color='red', lw=2)
import matplotlib.pyplot as plt
plt.figure(figsize=(15,15))
plt.scatter(phot_table4['residual_aperture_sum'], phot_table4['aperture_sum_0'])
plt.xlabel("Circular Apertures")
plt.ylabel("Background")
plt.xlim(100000,200000)
plt.ylim(150000,240000)
(150000, 240000)
pos = Table(names=['x_0', 'y_0'], data=[sources['xcentroid'],
sources['ycentroid']])
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from photutils.datasets import (make_random_gaussians_table,
make_noise_image,
make_gaussian_sources_image)
from photutils.psf import (IterativelySubtractedPSFPhotometry,
BasicPSFPhotometry)
from photutils import MMMBackground
from photutils.psf import IntegratedGaussianPRF, DAOGroup
from photutils.detection import DAOStarFinder
from photutils.detection import IRAFStarFinder
from astropy.table import Table
from astropy.modeling.fitting import LevMarLSQFitter
from photutils.psf import IntegratedGaussianPRF, DAOGroup
from photutils.background import MMMBackground, MADStdBackgroundRMS
from photutils.psf import IterativelySubtractedPSFPhotometry
from astropy.stats import gaussian_sigma_to_fwhm
sigma_psf = 2.0
bkgrms = MADStdBackgroundRMS()
std = bkgrms(data)
iraffind = IRAFStarFinder(threshold=3.5*std,
fwhm=sigma_psf*gaussian_sigma_to_fwhm,
minsep_fwhm=0.01, roundhi=5.0, roundlo=-5.0,
sharplo=0.0, sharphi=2.0)
daogroup = DAOGroup(2.0*sigma_psf*gaussian_sigma_to_fwhm)
mmm_bkg = MMMBackground()
fitter = LevMarLSQFitter()
psf_model = IntegratedGaussianPRF(sigma=sigma_psf)
photometry = BasicPSFPhotometry(group_maker=daogroup,
bkg_estimator=mmm_bkg,
psf_model=psf_model,
fitter=LevMarLSQFitter(),
fitshape=(11,11))
result_tab = photometry(image=data, init_guesses=pos)
result_tab
WARNING: The fit may be unsuccessful; check fit_info['message'] for more information. [astropy.modeling.fitting]
| x_0 | y_0 | flux_0 | id | group_id | x_fit | y_fit | flux_fit | flux_unc | x_0_unc | y_0_unc |
|---|---|---|---|---|---|---|---|---|---|---|
| float64 | float64 | float64 | int64 | int64 | float64 | float64 | float64 | float64 | float64 | float64 |
| 0.91894685 | 840.89515 | 104622.98672550803 | 1 | 1 | -0.3173792241660808 | 840.8844099643336 | 228890.43057933258 | 11388.231506821816 | 0.14968072965209928 | 0.09072328854096326 |
| 1.3110953 | 603.52503 | 261117.40290235524 | 2 | 2 | 1.300569260807128 | 603.1321144660371 | 320870.1218047408 | 9054.56657206123 | 0.09383415788255292 | 0.0787160379628428 |
| 1.3819622 | 194.14896 | 237618.21822195407 | 3 | 3 | 1.3695821448723215 | 194.1425495250382 | 303854.064838115 | 5582.723739886595 | 0.06073066600620981 | 0.05136082396918347 |
| 13.709729 | 966.33218 | 117309.86477837292 | 4 | 4 | 13.712264141803514 | 966.2986230219999 | 159191.53701778722 | 6206.966215302538 | 0.11151032905848834 | 0.11157241528485883 |
| 14.85393 | 104.58858 | 93019.76138578902 | 5 | 5 | 14.891850631606813 | 104.72455897028448 | 131748.69588285204 | 6487.774791122306 | 0.14084037569064306 | 0.14086124526028662 |
| 15.308049 | 662.25654 | 162433.68401351475 | 6 | 6 | 16.56388179248046 | 661.5843508658699 | 184976.22059674095 | 14269.97078440453 | 0.22297797670188038 | 0.220723213395023 |
| 26.026207 | 990.82951 | 89210.75882095125 | 7 | 7 | 26.02352000642823 | 990.787169671028 | 126253.19713086892 | 6245.721475319197 | 0.14145305270483385 | 0.14152292535230906 |
| 31.916545 | 171.79805 | 100658.13680020157 | 8 | 8 | 31.917581972682235 | 171.8343283450299 | 141613.77263797342 | 6448.423787265623 | 0.13023237686673134 | 0.13023927840828783 |
| 35.353299 | 267.10879 | 103527.16275763283 | 9 | 9 | 35.350323440435055 | 267.06805584474233 | 143578.56457971287 | 6366.09245212229 | 0.12684491251918442 | 0.12681320427197512 |
| 37.115099 | 83.290352 | 69293.6427486517 | 10 | 10 | 37.106416105574795 | 83.2152158858961 | 92957.74036971448 | 5166.642302788163 | 0.15896480832543164 | 0.15897906232912956 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 974.21449 | 695.40998 | 83716.48918003912 | 199 | 195 | 974.1577370810207 | 695.4605601376018 | 119046.20068618852 | 6389.634909089663 | 0.15353945315930215 | 0.15358879510160484 |
| 981.81682 | 32.809747 | 64075.48566281513 | 200 | 196 | 981.8147584053027 | 32.80930777250456 | 94197.52428649586 | 5928.216869567888 | 0.18003712687274853 | 0.17996729945719703 |
| 1014.2458 | 456.11819 | 115649.23680101393 | 201 | 197 | 1014.2525036848132 | 456.0948214803198 | 151671.88055892254 | 6114.457602676671 | 0.11534710689773145 | 0.1152792350735681 |
| 1023.7488 | 564.6262 | 85052.00460631734 | 202 | 198 | 1023.7318537578201 | 564.7185501263821 | 117954.2302764824 | 5902.815571003013 | 0.14316200088663744 | 0.14312117353528225 |
| 1039.4274 | 492.84468 | 107735.28446978563 | 203 | 199 | 1039.482610771519 | 492.79794173036663 | 143473.7446237493 | 6066.073862686504 | 0.12103845104603224 | 0.12090424187403902 |
| 1041.6488 | 101.70837 | 86397.26413256068 | 204 | 200 | 1041.6762768184842 | 101.69553365408699 | 119421.10767637611 | 6064.947953854878 | 0.1453071072326163 | 0.14524493521378323 |
| 1043.7634 | 119.08471 | 71455.85762679712 | 205 | 201 | 1043.7819088676713 | 119.1137619317121 | 103288.8871070871 | 6018.555237232858 | 0.16669977440373449 | 0.16661564847624116 |
| 1044.666 | 632.90166 | 77439.81192077725 | 206 | 202 | 1044.7029909726716 | 632.9280085526168 | 102270.05699210084 | 5505.366866039239 | 0.15400224655084926 | 0.153935963968129 |
| 1044.7256 | 804.82416 | 72539.8043643901 | 207 | 203 | 1044.728165492963 | 804.8553946343857 | 98589.33826088342 | 5259.265893653711 | 0.15259818086073976 | 0.1525560134208171 |
| 1048.2053 | 380.56477 | 114486.67845913697 | 208 | 204 | 1048.1557223219777 | 380.56131251673327 | 149602.08111747573 | 5893.051612545279 | 0.1126984174351405 | 0.11268198008765348 |
#norm = ImageNormalize(strecth=SqrtStrecth())
plt.figure(figsize=(15,15))
plt.scatter(phot_table4['aperture_sum_0'], result_tab['flux_fit'])
plt.xlabel("Aperture Sum")
plt.ylabel("PSF Value")
plt.xlim(150000,240000)
#plt.ylim(150000,240000)
(150000, 240000)
plt.figure(figsize=(30,30))
plt.subplot(1, 2, 1)
plt.imshow(data, cmap='viridis', aspect=1,
interpolation='nearest', origin='lower')
plt.title('Simulated data')
plt.colorbar(orientation='horizontal', fraction=0.046, pad=0.04)
plt.subplot(1 ,2, 2)
plt.imshow(residual_image, cmap='viridis', aspect=1,
interpolation='nearest', origin='lower')
plt.title('Residual Image')
plt.colorbar(orientation='horizontal', fraction=0.046, pad=0.04)
<matplotlib.colorbar.Colorbar at 0xb1b691128>